PDEs and Finite Elements (Hands on)

Oliver Rübenkönig

Wolfram Research Numerical Computation Development Team

Documentation

Guide page

Finite Element Method Guide

Tutorials

Solving Partial Differential Equations with Finite Elements

Element Mesh Generation

Element Mesh Visualization

Finite Element Programming

Reference pages

DirichletCondition — Dirichlet conditions
NeumannValue — NeumannValues
ToBoundaryMesh — convert various input to a boundary mesh
ToElementMesh — convert various input to a full mesh

Why finite elements

Arbitrary shaped domains

Large class of PDEs

The coefficient form

Consider the single linear partial differential equation

m((∂^2)/(∂(t^2)))u+d(∂/(∂t))u+∇·(-c ∇u -α u + γ ) +β·∇u + a u = f on  Ω

Specific PDEs

Coefficient Form: m((∂^2)/(∂(t^2)))u+d(∂/(∂t))u+∇·(-c ∇u -α u + γ ) +β·∇u + a u = f

Laplace:                                                                       ∇·(-c ∇u -α u + γ) +β·∇u + a u = 0

Poisson :                                                                     ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Helmholtz:                                                                 ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Convection/Diffusion/Reaction:      ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Heat:                                        ((∂^2)/(∂(t^2)))u+(∂/(∂t))u+ ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Wave:                                     ((∂^2)/(∂(t^2)))u+(∂/(∂t))u+ ∇·(-c ∇u -α u + γ) +β·∇u + a u = 0

What is needed to solve a PDE with the FEM

PDE - the physical setting

Region - the geometry

Boundary Conditions - the connection between PDE and everything outside of geometry

Poisson equation

-∇·(∇u)= f

PDE

pde=-D[u[x,y],{x,2}]-D[u[x,y],{y,2}]==1

Region

Ω=Rectangle[{0,0},{1,1}];

RegionPlot[True,{x,y}∈Ω]

Boundary conditions

(Subscript[Γ, D])={DirichletCondition[u[x,y]==0,True]};

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈Ω]

Plot

Plot3D[uif[x,y],{x,y}∈Ω]

Poisson equation - Exercise

1 Change the Region (e.g. a Disk)

Solution

Region

Ω=Disk[];

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈Ω]

Plot

Plot3D[uif[x,y],{x,y}∈Ω]

2 Change the PDE (e.g. a Helmholtz equation)

Solution

-∇·(∇u)+u= f

PDE

pde=-D[u[x,y],{x,2}]-D[u[x,y],{y,2}]+100u[x,y]==1

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈Ω]

Plot

Plot3D[uif[x,y],{x,y}∈Ω]

3 Change the Boundary Condition to be valid on part of the boundary

Solution

Boundary conditions

(Subscript[Γ, D])={DirichletCondition[u[x,y]==0,x≥0]};

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈Ω]

Plot

Plot3D[uif[x,y],{x,y}∈Ω]

Summary

Choose a PDE

Choose a Region

Specify boundary conditions

DirichletConditions specify the value of the dependent variable at the boundary

Neumann boundary values

Helmoltz:      ∇·(-c ∇u -α u + γ ) +β·∇u + a u = f

Region

Ω=Disk[];

PDE operator

op=-D[u[x,y],{x,2}]-D[u[x,y],{y,2}]+100u[x,y]-1

Boundary conditions

(Subscript[Γ, D])={DirichletCondition[u[x,y]==0,x≥0]};

Solve

uif=NDSolveValue[{op==0,(Subscript[Γ, D])},u,{x,y}∈Ω]

This makes a plot of the solution.

Plot3D[uif[x,y],{x,y}∈Ω]

Neumann boundary values - theory

The finite element method is based on the weak form of the differential equation. This form is obtained by taking the

Stationary Coefficient Form:∇·(-c ∇u -α u + γ ) +β·∇u + a u = f

multiplying it by a so-called test function ϕ, and integrating over the region Ω:

(∫_Ω)∇·(-c ∇u-α u+γ) ϕ+β·∇u ϕ+a u ϕ-f ϕdx=0 .

Integration by parts gives:

(∫_Ω)(c ∇u+α u-γ)·∇ϕ+β·∇u ϕ+a u ϕ-f ϕ dx-(∫_(∂Ω))(OverVector[n])·(c ∇u+α u-γ)ϕdx=0 .

This process is done internally. The integrand in the boundary integral (OverVector[n])·(c ∇u-γ+α u) can be replaced with g-q u and yields:

(∫_Ω)(c ∇u+α u-γ)·∇ϕ+β·∇u ϕ+ a u ϕ-f ϕ dx-(∫_(∂Ω))(g-q u)ϕdx=0 .

g-q u specifies the value of the boundary integral integrand (OverVector[n])·(c ∇u-γ+α u) in the weak form and thus the name NeumannValue. Other boundary conditions are conceivable, but currently not implemented.

NeumannValue boundary value - natural condition

(OverVector[n])·(c ∇u + α u - γ)= g - q u  on  ∂Ω

(Subscript[Γ, N])=NeumannValue[0,x<0]

This solves the equation:

uif=NDSolveValue[{op==(Subscript[Γ, N]),(Subscript[Γ, D])}, u,{x,y}∈Ω]

This makes a plot of the solution.

Plot3D[uif[x,y],{x,y}∈Ω]

NeumannValue boundary value - pure Neumann condition

(OverVector[n])·(c ∇u + α u - γ)= g - q u  on  ∂Ω

(Subscript[Γ, N])=NeumannValue[-1,x<0]

This solves and visualizes the equation:

uif=NDSolveValue[{op==(Subscript[Γ, N]),(Subscript[Γ, D])}, u,{x,y}∈Ω];
Plot3D[uif[x,y],{x,y}∈Ω,PlotRange->All]

NeumannValue boundary value - generalized Neumann condition - Robin condition

(OverVector[n])·(c ∇u + α u - γ)= g - q u  on  ∂Ω

(Subscript[Γ, N])=NeumannValue[2u[x,y],x<0]

This solves and visualizes the equation:

uif=NDSolveValue[{op==(Subscript[Γ, N]),(Subscript[Γ, D])}, u,{x,y}∈Ω];
Plot3D[uif[x,y],{x,y}∈Ω,PlotRange->All]

NeumannValue - Exercise

1 Axis of symmetry

Specify a symmetry region (e.g. a rectangular region with a hole: ).

Specify a Laplace equation.

Set multiple Dirichlet conditions u(x,y)==100  on the inner and u(x,y)==20  on the outer boundary.

Find a symmetric subsection  and compute the solution over that region. What are your findings? Why is this a big deal?

Solution

Neumann boundary conditions can be used to exploit symmetries of a region.

Specify a region with a hole:

Ω1=ImplicitRegion[-1≤x≤1&&-1≤y≤1&&!(-1/2<x<1/2&&-1/2<y<1/2),{x,y}];
RegionPlot[Ω1]

Specify a Laplace operator:

op=(Subscript[∇, {x,y}]^2)u[x,y];

Set multiple Dirichlet conditions u(x,y)==100  on the inner and u(x,y)==20  on the outer boundary:

Specify boundary conditions:

bcs={
DirichletCondition[u[x,y]==100.,Abs[x]==1/2&&-1/2≤y≤1/2||-1/2≤x≤1/2&&Abs[y]==1/2],
DirichletCondition[u[x,y]==20.,Abs[x]==1||Abs[y]==1]};

Solve the equation over a region and visualize the result:

uif1=NDSolveValue[{op==0,bcs},u,{x,y}∈Ω1];
Plot3D[uif1[x,y],{x,y}∈Ω1,PlotRange->All]

Specify a symmetric sub region:

Ω2=ImplicitRegion[x≤1&&x≥y&&y≥0&&(2 x≥1||2 y≥1),{x,y}];
RegionPlot[Ω2,PlotRange->{{-1,1},{-1,1}}]

Zero Neumann conditions are assumed for boundaries with no explicit condition set:

Solve the equation over a region and visualize the result:

uif2=NDSolveValue[{op==0,bcs},u,{x,y}∈Ω2];
Plot3D[uif2[x,y],{x,y}∈Ω2,PlotRange->{{-1,1},{-1,1}}]

Neumann boundary conditions can be used to exploit symmetries of a region.

Why is this a big deal? - Because the element / region of interest ratio is much higher, which may lead to a more accurate result.

ToElementMesh[Ω1]

ToElementMesh[Ω2]

2 Neumann flux over an edge

Here is a region specification

We would like to compute the temperature distribution in a pipe with cooling ducts:

Specify a region by boundaries:

boundaries={-y,(1/25)-(((-(3/2)+x))^2)-(y^2),1-(x^2)-(y^2), (x^2)+(y^2)-4,y-x*Tan[π/8]};
Ω=ImplicitRegion[And@@(#≤0&/@boundaries),{x,y}];

Show the region and the boundaries:

Show[RegionPlot[Ω],ContourPlot[Evaluate[Thread[boundaries==0]],{x,0.8,2.2},{y,-0.2,1.},ContourStyle->{Purple,Green,Red,Blue,Purple}],PlotRange->{{0.8,2.2},{-0.2,1.}},AspectRatio->Automatic]

Show the computational region in the overall region:

Show[Graphics[{Circle[{0,0},2],Circle[{0,0},1],Circle[#,1/5]&/@{{-1,-1},{-1,1},{1,-1},{1,1},{3/2,0},{-3/2,0},{0,3/2},{0,-3/2}}}],RegionPlot[Ω]]

Solve the Laplace equation with a thermal conductivity of 10 W/(m K) the temperature set to (200^o)C inside the pipe and (15^o)C out side the pipe.

Next, model a cooling liquid in the duct. The liquid absorbs -1000 (X?).

Solution

Solve the equation with only DirichletConditions.

uif=NDSolveValue[{-10 (Subscript[∇, {x,y}]^2)u[x,y]==0,{DirichletCondition[u[x,y]==200.,boundaries[[3]]==0.],
DirichletCondition[u[x,y]==15.,boundaries[[4]]==0.]}},u,{x,y}∈Ω,Method->{"FiniteElement","MeshOptions"->{"BoundaryMeshGenerator"->"Continuation"}}];

The "Continuation" mesh generator is used to better resolve the corners and edges of the domain.

Visualize the solution.

ContourPlot[uif[x,y],{x,y}∈Ω,Mesh->None,ColorFunction-> "TemperatureMap",Contours->21,AspectRatio->Automatic,PlotRange->All]

Solve a Laplace equation with a temperature set on the red and blue boundaries and a flux over the green boundary:

Solve the equation with a liquid in the cooling duct.

uif=NDSolveValue[{-10 (Subscript[∇, {x,y}]^2)u[x,y]==NeumannValue[-1000.,boundaries[[2]]==0.],{DirichletCondition[u[x,y]==200.,boundaries[[3]]==0.],
DirichletCondition[u[x,y]==15.,boundaries[[4]]==0.]}},u,{x,y}∈Ω,Method->{"FiniteElement","MeshOptions"->{"BoundaryMeshGenerator"->"Continuation"}}];

Visualize the solution.

ContourPlot[uif[x,y],{x,y}∈Ω,Mesh->None,ColorFunction-> "TemperatureMap",Contours->21,AspectRatio->Automatic]

Summary

NeumannValues can be used to model axis of symmetry

NeumannValues can be used to model flux over an edge (2D) or surface (3D)

Transient equations

Setup a region

Ω=ImplicitRegion[!((x-5)^2+(y-5)^2≤3^2),{{x,0,5},{y,0,10}}];
RegionPlot[Ω,AspectRatio->Automatic]

The dependent variable now depends on time and space.

Set up the PDE operator, being sure to include the time t as an argument of u.

op=(∂_t)u[t,x,y]-(InactiveLaplacian[u[t,x,y], {x,y}, Inactive])-20;

Also, the boundary conditions may depend on space and time.

Set up the boundary conditions.

(Subscript[Γ, D])=DirichletCondition[u[t,x,y]==0,x==0&&8≤y≤10];

(Subscript[Γ, N])=NeumannValue[-2*u[t,x,y],(-5+x)^2+(-5+y)^2==9];

To solve first-order time-dependent problems, NDSolve needs an initial condition and a time region.

Solve the PDE.

uifHeat=NDSolveValue[{op==(Subscript[Γ, N]),(Subscript[Γ, D]),u[0,x,y]==0},u,{t,0,100},{x,y}∈Ω];

Make an animation of the solution.

Animate[Plot3D[uifHeat[t,x,y],{x,y}∈uifHeat["ElementMesh"],PlotRange->{0,160},ColorFunction->"TemperatureMap",Mesh->All],{{t,10},0,40,2},SaveDefinitions->True]

Transient equations - Exercise

1 Heat equation

Use a transient version of the equation from above and ramp up the amount of heat absorbed for 10 seconds from 0 to -1000 and run the simulation for 25 seconds for the   region.

As the initial condition for the transient PDE the stationary solution without a NeumannValue can be used.

Solution

Specify a region by boundaries:

boundaries={-y,(1/25)-(((-(3/2)+x))^2)-(y^2),1-(x^2)-(y^2), (x^2)+(y^2)-4,y-x*Tan[π/8]};
Ω=ImplicitRegion[And@@(#≤0&/@boundaries),{x,y}];

Create the initial conditions from the stationary solution.

uinit=NDSolveValue[{-10 (Subscript[∇, {x,y}]^2)u[x,y]==0,{DirichletCondition[u[x,y]==200.,boundaries[[3]]==0.],
DirichletCondition[u[x,y]==15.,boundaries[[4]]==0.]}},u,{x,y}∈Ω,Method->{"FiniteElement","MeshOptions"->{"BoundaryMeshGenerator"->"Continuation"}}];

Define transient DirichletCondition and NeumannValues.

ΓD[t_]={DirichletCondition[u[t,x,y]==200.,boundaries[[3]]==0.],DirichletCondition[u[t,x,y]==15.,boundaries[[4]]==0.]};
ΓN[t_]=NeumannValue[If[t≤10.,-100. t,-1000.],boundaries[[2]]==0.];

Solve the transient PDE.

uif=NDSolveValue[{D[u[t,x,y],t]-10 (Subscript[∇, {x,y}]^2)u[t,x,y]==ΓN[t],ΓD[0],u[0,x,y]==uinit[x,y]},u,{t,0,25},{x,y}∈uinit["ElementMesh"]];

Visualize the solution.

framesNTV=Table[ContourPlot[uif[t,x,y],{x,y}∈uif["ElementMesh"],Mesh->None,ColorFunction->"TemperatureMap",Contours->21,PlotRange->All,AspectRatio->Automatic,Frame->False]//Quiet,{t,1/2,15,1/2}];
Manipulate[framesNTV[[i]],{{i,9,"time"},1,Length[framesNTV],1},SaveDefinitions->True]

Interface - API - Documentation

All low level functions are documented!

Initialization
InitializePDECoefficients — initialize partial differential equation coefficients
InitializeBoundaryConditions — initialize boundary conditions
InitializePDEMethodData — initialize partial differential equation method data

PDECoefficientData ▪ BoundaryConditionData ▪ FEMMethodData

Discretization
DiscretizePDE — discretize initialized partial differential equations
DiscretizeBoundaryConditions — discretize initialized boundary conditions

DiscretizedPDEData ▪ DiscretizedBoundaryConditionData

Solution
DeployBoundaryConditions — deploy discretized boundary conditions into discretized partial differential equations
LinearSolve — solve linear systems of equations

Tutorial

Finite Element Programming

Guide page

Finite Element Method Guide

Thanks and - NDSolve - Rocks!

Summary:

Linear, variable coefficient, transient systems of PDEs over arbitrary regions in
1D/2D/3D

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;